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Abstract. Two local discontinuous Galerkin (LDG) methods using some non- 
standard numerical fluxes are developed for the Helmholtz equation with the 
first order absorbing boundary condition in the high frequency regime. It 
is shown that the proposed LDG methods are absolutely stable (hence well- 
posed) with respect to both the wave number and the mesh size. Optimal order 
(with respect to the mesh size) error estimates are proved for all wave numbers 
in the preasymptotic regime. To analyze the proposed LDG methods, they are 
recasted and treated as (non-conforming) mixed finite element methods. The 
crux of the analysis is to establish a generalized inf-sup condition, which holds 
without any mesh constraint, for each LDG method. The generalized inf-sup 
conditions then easily infer the desired absolute stability of the proposed LDG 
methods. In return, the stability results not only guarantee the well-posedness 
of the LDG methods but also play a crucial role in the derivation of the error 
estimates. Numerical experiments, which confirm the theoretical results and 
compare the proposed two LDG methods, are also presented in the paper. 



1. Introduction 



This paper is the third installment in a series (TOJ [IT] which devote to developing 
absolutely stable discontinuous Galerkin (DG) methods for the following prototyp- 
ical Helmholtz problem with large wave number: 

(1.1) -Au-k 2 u = f in il C K d , d = 2,3, 

Ou 

(1.2) Viku = g onr = dfl, 

on n 

where i = \] — 1 denotes the imaginary unit, k 6 R + is a given positive (large) 



number and known as the wave number. (1.2) is the so-called first order absorbing 
boundary condition [9]. 

We recall that |10[ lllj focused on designing and analyzing h- and hp-interior 
penalty discontinuous Galerkin (IPDG) methods which are absolutely stable (with 
respect to wave number k and mesh size h) and optimally convergent (with respect 
to h) . The main ideas of [101 lllj are to introduce some novel interior penalty terms 
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in the sesquilinear forms of the proposed IPDG methods and to use a non-standard 
analytical tool, which is based on a Rellich identity technique, to prove the desired 
stability and error estimates. The numerical experiment results shown that the ab- 
solutely stable IPDG methods significantly outperform the standard finite element 
and finite difference methods, which are known only to be stable under stringent 
mesh constraints hk < I or hk 2 < 1 (cf. 8, 16 ), for the Helmholtz problem. More- 
over, the numerical experiment results also shown that these IPDG methods are 
capable to correctly track the phases of the highly oscillatory waves even when the 
mesh violates the "rule-of-thumb" condition (i.e., 6—10 grid points must be used 
in a wave length). The main difficulty of analyzing the Helmholtz type problems is 
caused by the strong indefinitcness of the Helmholtz equation which in turn makes 
it hard to establish stability estimates for its numerical approximations. The loss of 
stability in the case of large wave numbers results in an additional pollution error 
(besides the interpolation error) in the global error bounds. Extensive research has 
been done to address the question whether it is possible to reduce the pollution 
effect, we refer the reader to Chapter 4 of [15] and the references therein for an 
detailed exposition in this direction. 

Motivated by the success of [TOl ITT] , the primary objective of this paper is to 
extend the work of [101 111] to the local discontinuous Galerkin (LDG) formulation, 
which is known to be more "physical" and flexible than the IPDG formulation on 
designing DG schemes [U [5] . As it is well-known now, the key step for constructing 
LDG methods is to design the numerical fluxes. As soon as the numerical fluxes are 
selected, for a large class of coercive elliptic and parabolic second order problems, 
there is a general framework for carrying out convergence analysis of LDG methods 
PQ. Unfortunately, this general framework does not apply to the Helmholtz type 
problems which is extremely noncoercive/indefinite for large wave number k. Nev- 
ertheless, when designing the numerical fluxes for our LDG methods, we borrow 
the idea of [1] by establishing the connection between our LDG methods and the 
IPDG methods of QH [TT] although it turns out that the IPDG methods of [TU1 E] 
do not have exactly equivalent LDG formulations due to the non-standard penalty 
terms used in [TUHU]. This then leads to the construction of our first LDG method. 
It is proved and numerically verified that this LDG method is absolutely stable and 
optimally convergent for the scalar variable. However, it is sub-optimal for the 
vector/flux variable. To improve the approximation accuracy for the vector/flux 
variable, we design another set of numerical fluxes which result in the construction 
of our second LDG method. It is proved that the second LDG method is also ab- 
solutely stable and gives a better approximation for the vector/flux variable than 
the first method. On the other hand, it is computationally more expensive than 
the first LDG method, which is expected. 

To analyze the proposed LDG methods, we take an opposite approach to that 
advocated in pQ, that is, instead of converting LDG methods to their "equivalent" 
IPDG methods in the primal form, we recast and treat our LDG methods as non- 
conforming mixed finite element methods. To avoid using the standard techniques 
such as Schatz argument (cf. [21 H]) or Babuska's inf-sup condition argument [16] 
to derive error estimates (and to prove stability), both approaches would certainly 
lead to stringent mesh constraints, our main idea is to establish a generalized inf-sup 
condition, which holds without any mesh constraint, for each LDG method. The 
generalized inf-sup conditions then immediately infer the desired absolute stability 
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of the proposed LDG methods. In return, the stability results not only guarantee 
the well-posedness of the LDG methods but also play a crucial role in the derivation 
of the (optimal) error estimates. 

It should be pointed out that a lot of work has recently been done on develop- 
ing DG methods using piecewise plane wave functions, oppose to simpler piecewise 
polynomial functions as done in this paper, for the Helmholtz type problems. How- 
ever, to the best of our knowledge, none of these plane wave DG method is proved 
to be absolutely stable with respect to wave number k and mesh size h. We refer 
the reader to [TH Q31 [T7] and the references therein for more discussions in this 
direction. We also refer to [TOl HJJ for more discussions and references on other 
discretization techniques for the Helmholtz type problems. 

This paper consists of four additional sections. In Section [2j we introduce the 
notations used in this paper and present the derivations of our two LDG methods. 
In Section [3j we present a detailed stability analysis for both LDG methods. The 
main task of the section is to prove a generalized inf-sup condition for each proposed 
LDG method. Similar to jTOl [IT], a nonorthodox test function trick is the key to 
get the job done. In Section [1J a non-standard two-step error estimate procedure 
is used to derive error estimates for the proposed LDG methods. Once again, the 
stability estimates established in Section [3] play a crucial role. Finally, Section [5] 
contains some numerical experiments which are designed to verify the theoretical 
error bounds proved in Section [4] and to compare the performance of the proposed 
two LDG methods. 

2. Formulation of local discontinuous Galerkin methods 

The standard space, norm and inner product notation are adopted in this paper. 
Their definitions can be found in [TJ [2 [H [JJ5] . In particular, (-,-)q and (-,-)s for 
£ C dQ denote the L 2 -inner product on complex-valued L 2 (Q) and L 2 (YA spaces, 
respectively. (•,•) := an d (•,•) := (- 7 -)dn- Throughout the paper, C is used 

to denote a generic positive constant which is independent of h and k. We also use 
the shorthand notation A < B and B > A for the inequality A < CB and B > CA. 
A ~ B is a shorthand notation for the statement A < B and B < A. 

Assume that Q C R d (d = 2, 3) is a bounded and strictly star- shaped domain with 
respect to a point xq G fi. We now recall the definition of star-shaped domains. 

Definition 2.1. Q C R d is said to be a star-shaped domain with respect to xq G Q 
if there exists a nonnegative constant cq such that 



Let 7ft be a family of partitions of fi parameterized by h > 0. For any tri- 
angle/tetrahedron K G 7ft, we define hx '■= diam(if) and h := maxx^Th dis- 
similarly, for each edge/face e of K G Th, define h e := diam(e). We assume that 
the elements of Th satisfy the minimal angle condition. We also define 



(2.1) 
Q c R d 



(x - Xq) ■ Uq > Cq Vx G dQ. 

is said to be strictly star-shaped if cq is positive. 
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Let e be an interior edge shared by two elements K\ and K 2 whose unit outward 
normal vectors are denoted by ri! and n 2 . For a scalar function v, let Vi = v\gKii 
and define 

W = 7^( v i + H =v K - v K >, [[«]] = uini + v 2 n 2 on e G ^ 

where K is i^i or if 2 ! whichever has the bigger global labeling and if' is the other. 
For a vector field v, let Vj = v\dKi and define 

W = ^( v i + v 2) ; [v]=vjf-vjf(, [[v]] = vi • ni, + v 2 • n 2 on e E £-{. 

As it is well-known now (cf. []]) that the first step for formulating an LDG 
method is to rewrite the given PDE as a first order system by introducing an 
auxiliary variable. For the Helmholtz problem (1.1)-(1.2) we have 

(2.2) (j = Vti in ft, 

(2.3) -diver ~k 2 u = f in ft, 

du 

(2.4) h lfcu = g on T, 

Clearly, the vector-valued function (often called the flux variable) er is the auxiliary 
variable. 

Then, multiplying (2.2) and (2.3| by test functions t and v, respectively, and 
integrating both equations over an element K E Th yields 



(2.5) / er • t dx = — / udivTdx+ / uxiK-Tds, 

Jk Jk JdK 

(2.6) / cr-\7vdx — k 2 / uv dx = fvdx+j cr-tiKvds, 
Jk Jk Jk JdK 

where \\k denotes the unit outward normal vector to dK. The above equations 
form the weak formulation one uses to define LDG methods for the Helmholtz 
problem (D}-jL2 ). 

Next, we define LDG spaces as follows 

V h := {v e L 2 (Q); Re(v)\ K ,Im(v) G P r (K) VK e %}, 

S ft := {t e (L 2 (n)) d ; Re(r)\ K , Im(r) e (Pf(^)) d VK € T h }, 

where P r (K) (r > 1) stands for the set of all polynomials of degree less than or 
equal to r on K . 

Finally, we are ready to define the following general LDG formulation: Find 
(uh,<Th) S Vh x S/j such that for all K E Th there hold 



(2.7) / a h ■ T h dx = u h divTft, dx + u K x\ K - T h ds, 

Jk Jk Jok 

(2.8) / cr h ■ Vv h dx-k 2 I u h v h dx = I fv h dx + a K ■ n K v h ds 
Jk Jk Jk Jok 

for any {vh,Th) € Vh x S^. Where the quantities and &k, which are called 
numerical fluxes, are respectively approximations to er = V« and u on the bound- 
ary dK of K. As it is well-known now that the most important issue for all LDG 
methods is how to choose the numerical fluxes. The different choices of the numer- 
ical fluxes obviously lead to different LDG methods. It is easy to understand that 
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these numerical fluxes must be chosen carefully in order to ensure the stability and 
accuracy of the resulted LDG methods. 

In this paper, we shall only consider the linear element case (i.e., r = £ = 1) and 
propose two sets of numerical fluxes (uk,o'k), which lead to two LDG methods. 
Our choices of numerical fluxes are inspired by the interior penalty discontinuous 
Galerkin (IPDG) methods proposed by Feng and Wu [10] and are identified with the 
help of the unified DG framework of [T] which bridges the primal DG formulations 
(e.g. IPDG methods) and the flux DG formulations (i.e., LDG methods). 

(1) LDG method #1: Set 

(T K = {V h u h } - i/3[[u h }}, u K = {v-h} + ^[[VhUh]] on e G £[, 
&k = -lfeifciiK + gtiK, uk = Uh on e G £;f . 

(2) LDG method #2: Set 

<trt = {<Th} - i/3[[tth]], u K = {u h } + \5[[a h ]] on e G £[, 
& K = -iku h n K + gn K , u K = u h on e G £,f . 

Where /3 and 5 are positive constants to be specified later and Vh denotes the 
piecewisely defined gradient operator over Th, that is, 'Vhlx = ~S/\k Vlf G Th- 

For the reader's convenience, we now briefly sketch the derivation of the primal 
DG formulation corresponding to our LDG method #1 by adapting the derivation 
given in the general framework of pQ. 

Substituting the numerical fluxes of LDG method #1 into (2.7 1 and (2.8), sum- 
ming the resulting equations over all element K G Th and using the following 
integration by parts identity 

{u h ,divT h ) n = -(VhUh,Th)n + ^ ((W}> HMDe 

+ ([ML {Th})e) + (uh,n K -T h ) Tl 

we get 

(2.9) (cr h , Th)n ~ {^hUh-r h )n 

, ... ) =0, 



( lS ([[VhUh]}, [[Th]]) e ~ <[[«/.]], {r h }) t 
(2.10) (<r h , V h v h ) n - fc 2 K, v h ) n - ^ ({V h u h } - i/3[[u h }}, [[v h ]]) e 



+ ifc {uh, Vh) r = (/> + (#= u /i)r • 



Setting t — Vu/s in (2.9) and subtracting the resulting equation from (2.10) then 



leads to the following formulation: 



(2.11) 



A h (u h ,v h ) - k 2 (u h ,v h ) n = F(v h ) Vv h G V h , 



6 



XIAOBING FENG AND YULONG XING 



where 



A h (u h ,v h ) 



(V 'hUh,V 'hV h )n + ik (u h ,v h ) r 
^ i(6([[V h u h ]],[[V h v h ]]) ( 



E 



^(N],W)e 



(([[«h]],{v h « h }) e + ({v h u h }, [K] 



^K) := (f,Vh)n + {g,v h ) l 



Hence, (2.11) is the corresponding primal (i.e., IPDG) formulation of our LDG 
method #1. Comparing (2.11) with the IPDG formulation of [TU], we see that all 
the interior penalty terms of (2.11 ) also appear in the IPDG formulation of [10] . In 
fact, it is exactly by reversing the order of the above derivation that leads to the 
discovery of the numerical fluxes of LDG method #1. 

3. Stability analysis 
From their constructions, it is easy to see that both LDG methods proposed in 



the previous section are consistent schemes for the Helmholtz problem (1.1 )— ( 1.2 ) 



For coercive elliptic (and parabolic) problems, the stability of such a numerical 
scheme can be proved easily as demonstrated in [1] (the same statement is true 
for their corresponding PDE stability analysis). However, the Helmholtz problem 



(1.1 )-( 1.2 ) is an indefinite problem and it becomes notoriously non-coercive for large 



wave number k. Deriving its stability estimates (i.e., a priori estimates of its PDE 
solution), particularly wave-number-dependent estimates, has been proved not to be 
an easy job (cf. [6j [7l [T3l HB] and the references therein). Numerically, such a quest 
has been known to be even harder because of the lower order of the smoothness 
and the inflexibility of (piecewise) approximation functions (cf. [5J [TFJ US] and 
the references therein). The stability of the numerical methods in all the above 
quoted references was proved under some very restrictive mesh constraints. An 
open question was then raised by Zienkiewicz [20] which asks whether it is possible 
to construct an absolutely stable (and optimally convergent) numerical method 
(i.e., no restriction on the mesh size h and the wave number k) for the Helmholtz 
equation. Almost a decade later Feng and Wu (TUJ E] were able to design for the 
first time such numerical methods, which happen to belong to the IPDG family, for 



the Helmholtz problem (1.1 )— ( 1.2 ) 



The goal of this section is to show that in the case of the linear element (i.e., 
r = t = 1) the LDG method #1 and #2 proposed in the previous section have 
the same stability property as the IPDG methods of [TU1 E] do, that is, the LDG 
method #1 and #2 are absolutely stable for all mesh size h > and all wave number 
k > without imposing any constraint on them. To establish this result, unlike the 
approach used and advocated in the general framework of ;.lj which converts an LDG 
method into its equivalent primal method and then analyzes the latter using the 
standard finite element Galerkin techniques, we shall directly fit the LDG method 
#1 and #2 into the (nonconforming) mixed method framework as done in |16| and 
adapting the mixed method techniques to prove the desired stability result. It is 
well-known that the key ingredient for the mixed method approach is to establish 
the inf-sup or Babuska-Brezzi condition for the (augmented) sesquilinear forms for 
each method. However, we are not able to prove such an inf-sup condition for cither 
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method without imposing mesh constraints (which we believe is not possible). To 
overcome the difficulty, our main idea is to prove a generalized (and weaker) inf-sup 
condition which holds for all h, k > 0. It turns out that this generalized inf-sup 
condition is sufficient for us to establish the desired absolute stability for the LDG 
method #1 and #2. To prove the generalized inf-sup condition, the key technique 
we use is a special test function technique which was first introduced and developed 
in [TO]. 

3.1. Absolute stability of LDG method #1. We first recast the LDG method 
#1 as the following nonconforming mixed method, which is easily obtained by 
adding (2.9) and ( 2.10 |. Find (u h ,cr h ) e 14 x S/, such that 

(3.1) A h (u h ,<T h ;vh,Th) = F(v h ,Th) Vf^n) e 14 x Ej,, 
where 

(3.2) A h (w h , Xh'i v h, Th) = (X/j, V h v h ) n - k 2 (w h , v h )n + ik (w h , v h ) r 

£ ({V h w h }-iP[[w h ]],[[v H ]]) e 

(16 ([[V h w h ]], [[T h ]]) e ([M, {r h }) e 

(3.3) F(v h , Th) ■= (/, v h ) u + (g, v h ) r . 

3.1.1. A generalized inf-sup condition. The goal of this subsection is to show 
that the sesquilinear form Ah defined in (3.2) satisfies a generalized inf-sup con- 
dition, which will play a vital role for us to establish the absolute stability of the 
LDG method #1 in the next subsection. 

Proposition 3.1. There exists an h- and k- independent constant c\ > such that 
for any {w h ,Xh) e V h X 

Re A h (w h ,Xhi v h,Th) 



(3.4) sup 



ImAh(w h ,Xh, v h,T h ) ci 
+ Sup n n > —\\Wh\\DG, 

(v h ^ h )ev h xs: h \\Vh\\DG 7l 



where 



, \ , 8 fk 2 + l 1 1 

(3.5) 7i := 1 + k + \ - + max — — h — + — 

V 8 eefiV Bh e h 2 e Bh\. 

1 

(3.6) \\w h \\DG ■= (k 2 \\w h \\ 2 L 2 m + k 2 \\w h \\ 2 L2{r) + c n ||V/ l ^|||,2 (r) + \w h \l :h ) * . 

\w h \i, h ■= ( \\ Vw h\\h(K)) 2 ■ 
Ken 

Proof. The main idea of the proof is that for a fixed (wh,Xh) Vh x Sj, we need to 
pick up two sets of special functions {vh,T~h) € Vh X E/j, which, as expected, must 
depend on (wh,Xh) £ Vh x E/,, such that both quotients in (3.4) can be bounded 
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from below by ||w/i||dg- When that is done, the inf-sup constant ci/71 will be 
revealed in the process. Since the proof is very long, we divide it into four steps. 

Step 1: Taking the first test function. 

We first choose the test function (vh,Th) — (w/i, — V/jUi/j) to get 

Ah(Wh,Xhi v h,Th) = A h {w h ,Xh'i w h,~^hW h ) 

= (VhWh, S7hWh)n - k 2 (w h ,w h ) n + ik (w h ,w h ) r 

+ e (iWhWh]}, i[v h wh]}) e - ([K]], {Vhw h }), 



E ({Vhw h }-il3[[w h ]],[[w h ]]) t 



Taking the real and imaginary parts yields 
(3.7) ReA h (w h ,Xhl w h,-'VhWh) 

= \w h \l h - k 2 \\w h \\l 2(n) - 2 Re E (M, {VhW h }) t 



(3.8) Im A h (w h ,Xh^ w h,-^hWh) 



Step 2: Taking the second test function. 

Inspired by the special test function technique of [TU], we now choose another 
test function (vh, T/,) = (a ■ V/iiu/t, — V^w^) with a :— x — xq (see Definition 2.1 ) 
and use the fact that V h^h = ^h w h to get 

Ah(wh,Xh'i v h,i~h) = A h (w hl x h \ a ■ ^hWh,-^hWh) 

= {V h w h ,V h Wh)n - k 2 (w h ,a ■ V h w h ) n + ik (w h ,a ■ V h w h ) r 

+ E { l6 ([\?hw h ]], [[V h w h ]]) e - ([[w h ]},{V h w h }) ( 



E ({^hw h } - i(3[[w h ]}, [[a ■ V h w h ] 



ee£! 



Taking the real part immediately gives (note that Vh — ol ■ V ' hWh) 
(3.9) ReA h (w h ,Xh'- v h,-^hW h ) 

= \ w h\l,h ~ k 2 (w h ,v h ) n - kIm(w h ,v h } T - E ([ML [K]]) c 



E Rc ( ([[w h ]],{V h w h }) e + ({V h w h }, [[v h ]]), 



ees! 



Step 3: Deriving an upper bound for k 2 \\ WhW 2 ^^ ■ 

To bound (3.7) from below, we need to get an upper bound for the term k ||i0/i||r,2 
on the right hand side of (3.7). This will be done by carefully and judicially com- 
bining (3.8) and (3.9) with some other differential identities, which we now explain. 



(O) 
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Using the integral identity 



(3.10) 2\\w h \\ L 2( K) = / a-n K \w h \ ds - 2 Re(w h , v h ) K - (d - 2)\\w h \\ L 2 {K) , 

J dK 



and (3.7), (3.9 1 we get 



(3.11) 2k 2 \\w h \\ 2 L2{n) = 2fte A h (w h ,x h ,v h ,-Vhw h ) 
+ (d-2)EeA h (w h ,Xhi w h,-'^hWh) 
+ 2klm (w h ,v h ) r - 2\w h \l th - (d-2)\w h \l h 

+ 2 Re ((IK]]> {VhW h }) e + ({Vw h }, [K]] 
+ 2(d-2)Rc ([K]],{V h «> h }) e 



V(2^Im([K]],[M) e + fc 2 V / a-n K \w h \ 2 ds 
v i? c t Jbk 



By the elementary identity \a\ 2 — \b\ 2 = Re(a + b)(a — b) for any two complex 
numbers a and b we have 

(3.12) ^2 / a-n^lw/jpds = 2Re (a • n e {w h }, [w h ]) e + (a ■ n n , M 2 ) r • 



Next, using the local Rcllich identity (see [TOJ Lemma 4.1]) 

(d-2)\\Vw h \\ 2 LHK) +2Re(Vw h ,Vv h ) K = / a ■ n K \Wv h \ 2 ds 

JdK 

and the fact that Vh v h — V '/,%, we get 

(3.13) dK|^ = (d-2)K|^ + 2Re ^ (V«J h , Vv a )k 



ffef h 



Re ^ / a • n^Vv/il 2 ds 

KGTh ■' dK 



2 Re ^ (a • n e {V h u;/,}, [Vfcitf/,]) e + ^ (a • n e , |V /l w h | 2 ) £ 
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Substituting (3.12) and (3.131 into (3.11) we get 



(3.14) 2k 2 \\w h \\ 2 L 2 in) = 2Re A h {w h ,x h ;vh,-Vhw h ) 
+ {d-2)ReA h (w h ,x h '-, w h,-VhWh) 
+ 2k 2 Re (« ' ne{w h }, [w h ]) e + k 2 (a ■ n n , \w h \ 2 ) J 



ees! 



2fclm (w h ,v h ) r - ^ {a ■ n ei \\' h w h \ 2 ) e 
2 Re ((" ' n e {V h w h ], [V h w h ]) e - ({V h w h }, [[v h ]]) t 
2(d - 1) Re ]T ([[v>h]], {V h w h }) e + 21m £ /? ([[w h ]], [M> £ 



To get an upper bound for k W^hW^^iy we need to bound the terms on 
right-hand side of (3.14), which we now bound as follows. 



the 



(3.15) 2fc 2 Re (a-n e {w h },[w h }) e <Ck 2 £ h e h ||^|U=(K e u^)||K]|U=(e) 

eesi eeel 

< yiKiii^ + cE ^-mwH]\\h (e) 

(3.16) k 2 (a-n n ,\w h \ 2 ) r <Ck 2 \\ 

It follows from the star-shaped assumption on f2 that 

(3.17) 2klm{w h ,v h ) r - ^ (a • n e , \V h w h \ 2 ) e 



^ V(e) 



< Ck 2J ||w/ l ||L2 (e )||V/ l Wh|| L 2 (e ) - co 2J HV^Whl' 2 



cn 
2 

By the trace inequality [2], we also have 
(3.18) 2dReJ2{[[w h ]],{V h w h }) e 



<2d^2h e 2 ^2 \\^hW h \\ L 2 {K )\\[w h }\\ L 2 ie) 

e€ £i h K=K e ,K' e 

<lKI?, h + ^E ^IIWIIW)- 
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(3.19) - 2 Re E (<« • M^hWh}, [V h w h ]) e - ({V h w h }, [[v h ]]) e ) 



ee£! 



2Re^ 

eee' 



3 = 1 



E / ((" • HHVhWh ■ n e } 



-(a • n e ){V h w h ■ T 3 SjV h [wh] 



d-l 



<EE^ 5 E HVfcWfcii^wiiiVfcWfc-Tiiiu. 



(e) 



d-l 



< ^KlU + C E ^7 E ^III V ^ • ^elll W 



eE£j 



3 = 1 



By the definition of Vh := a • V^w/j, we get 



(3.20) 21m E /?<[K]MM>e = 21m E ^(WJ%])e 

/ d-l 

= 21m E Z 3 ( M> ( a ' n e )V ft w h • n e + E( a ' T 3 e )V h w h ■ t{ 

eeSl \ 3 = 1 

< C E ^11 Nfo]IU 2 (e) II [VftiOh. • n e ]|| L 2 (e) 

d-l 

+ C E ^11 KHU'(e) E II Fhw h ■ ri\\\ L , (e) 

eeS{ 3 = 1 



^ C \l f E (^11 WH^(e) +^II[V^ • ne]||| 2(e) ) 



C/3 E (ll M WUe) + E II l V ^ ' ri]\\h(e)) > 



where {t^}^ = ^ denotes an orthogonal tangential frame on the edge/face e, and we 
have used the decomposition a = (a ■ n e )n e + J2j=i( a ' T i) T i- 
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Now substituting estimates ( 3.15 )— ( 3.19 1 into (3.14) we obtain 

(3.21) 2k 2 \\w h \\l 2{n) <2~Rf}A h (w h ,Xh,Vh,-VhWh) 

+ (d - 2)ReA h (w h ,Xh'i w hi-^hWh) 



+ Ck 2 \\w h \\ 2 L 2 {T) - y ||Vh«; h ||| 3(r) + y||^|||2 (fi) 



C E ^IIK]H^(e)- 2Re E (M.{VfcWh}) e 
eesl P e ee£! 



+ C E V £ (^ Kl lli^e) + <*ll [Vfctflh ■ n e ]||| 2(e) 
+ C^(/?||K]||| 2(e) +^/3||[V^ 



d-1 

• T i]lll 2 (e) 

ee£! ' 3=1 



On noting that (3.8 ) provides upper bounds for terms ||[VfeWft-n e ]||^ 2 ( e \, || [w/J ||^2( e ) 
and fc 2 ||w^||^ 2( . r x in terms of Im-A/^u;/,, \ h ; Wh, — ^hWh), using these bounds in 



(3.21) we get 



(3.22) 2k 2 \\w h \\ 2 L2(n) + y|K||£ 2(r) + C f\\y h w h \\ 2 L2(r) 



< 2ReA h (w h ,Xhi v h,^hW h ) + (d- 2)ReA h (wh,Xhl w h,-^hWh) 



lk 2 . 

;Kli,fc + -d|wh||i» ( n) " 2Re E ({ V ^}= [KID. 



cE^ + ^E^ii^-^iiW) 



d-1 



e 6 £j[ 3=1 

.Mi Im^(iOh, Xh'^ w h, -V h w h ), 



where 



(3.23) Mi :=C(l + fc + W^+max^-^ 



To bound the jumps of the tangential derivatives HIV^u^ • Te]|lz,ss( e ) in (3.21), 
we appeal to the inverse inequality 

(3-24) \\[V h w h ■ ri]\\l He) < Ch- 2 \\[w h ]\\ 2 L2{e) , 
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and using (3.8) and (3.7) to get 



< 2Re A h (w h ,Xh\ v h, -V h^h) + {d ~ 2)ReA h (w h ,Xh-, w h,-^hWh) 

k 2 

+ Y W w h\\h(n) + M 2 lmA h (w h ,x h ;wh,-^hWh) 
~ 2 Re ({VhW h }, [K]]) e + \\w h \l h 

< 2ReA h (w h ,Xh': v h,-'VhW h ) + (d- 1) Re A h (w h , Xh'i w h , ~V h w h ) 

3fc 2 1 
+ + M 2 lmA h (w h ,Xh;w h ,~V h w h ) - -\w h \ 2 lh , 



where 



(3.25) ^ = ^ + -g(| + || 



nfi i If 3 /k 2 + l 1 1 \ 



Using the linearity of the sesquilinear form A^ we have 

(3.26) fc 2 |K|||2 (n) + A; 2 ||^||| 2(r) + c n \\V h w h \\ 2 L2{T) + \w h \l h 

< AReA h (w h ,x h ;vh; -V h w h ) +2(d- 1) Re A h (w h , Xhl w h, -V/.W/,) 

+ 2A^ 2 Im4K, Xh> w h,~^hW h ) 
= ReA h (wh,Xh'iU>h, -VhW h ) + 2M 2 lmA h (w hl Xh> w h, ~V h w h ). 

with Wh = 4u/i + 2(d — l)wh- 

Step 4- Finishing up. 



By the definition of || • \\dg m (3.6) and the fact that V^f/i = ^hWh we have 

(3.27) |K||1, G = fc 2 ||c* • V fe ^||| 2(n) + k 2 \\a ■ V h w h \\ 2 L2(T) 

+ c n \\\7 h w h \\ 2 L 2 (r) + \w h \\ h 
< Ck 2 \\V h w h \\ 2 LHn) + Ck 2 \\V h w h \\ 2 L2[v) 

+ c { -i\\V h w h \\ 2 L2{T) + \w h \\ h 
<C(l + k 2 ){\w h \l h + c n \\V h w h \\l 2(r) ] 
<C(l + k 2 )\\w h \\ 2 DG . 

Thus, it follows from the triangle inequality that 



(3.28) WwhWdg < 4K||dg + 2(d - l)\\wh\\DG < C(l + k)\\w h \\ DG . 
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Now from ( |3.26[ ) and (3.28) we have 

(3 29) ^ eA h ( w hiXh>^h,-Vh'Wh) , Im4K,X/ l ;«'/.>-V/,w/ i ) 



> 



> 



\\Wh\\DG \\Wh\\DG 

ReA h (w hl Xh'^h, -Vhw h ) lmA h (w h , Xh\ w h, ~VhW h ) 
C(l + k)\\w h \\ DG \\Wh\\DG 
1 ReA h (w h , Xh> Wh, + 2M 2 Im A h (w h , w h ,-V h w h ) 



2M 2 

> —\\wh\\DG 

7i 



\Wh\\DG 



for some constant ci > and ji is defined by (3.5). Hence, (3.4) holds. The proof 
is complete. □ 

Remark 3.1. (a) We note that 71 depends on both h and k. 

(b) The generalized inf-sup condition is a weak estimate because it does not pro- 
vide a control for the variable Xk- As a comparison, we recall that the standard 
inf-sup condition for the sesquilinear form Ah should be 



sup 

(v h ,T h )eV h y.H h 



\Ah(w h ,x h ;v h ,T h )\ 



>ci\\{w h ,x h )\\ V(w h ,Xh) G V h x H h 



\\{Vh,Th)\\ 

for some positive constant c\ = ci(k, (3, 5, fl). Where 

\\(wh,Xh)\\ ■= (k 2 \\w h \\ 2 L 2 (n) + WXhWhin)) 1 " ■ 

However, the above standard inf-sup condition can be proved only under the mesh 
constraint h — 0(k~ 2 ) and we believe that it does not hold without a mesh con- 
straint. 

3.1.2. Stability estimates. The goal of this subsection is to establish the absolute 
stability for the LDG method #1 using the generalized inf-sup condition proved in 
the previous subsection. 



Theorem 3.1. Let (u h ,er h ) e V h x T, h solve (jO). Define 

(3-30) M(/, ff ):=||/||n + NU»(r). 

Then there hold the following stability estimates: 

(3.31) K||dg <lik^M{f,g). 

(3.32) |KIU= ( fi) < lik-Hl + (S + r'Kmax h K 1 )) M(f,g). 
Proof. By Schwarz inequality we have 

(3.33) \F(v h ,x h )\ < ll/IU 2 (o)II^IU 2 (o) + HffllL 2 (r)IMlL 2 (r) 

< Ck-Hlif^g) (fc 2 ||^||i 2(0) + fc 2 |k||! 2(r) ) ; 



<Cfc- 1 M(/, ff )||^|| 



DG- 
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Let {wh,Xh) — { u h,Ch) in (3.4). By equation (3.1) and (3.33) we get 



^1 
7i 



\Uh\\DG < 



sup 



ReA h (u h ,a h ;v h ,T h ) 



sup 



sup 

(v h ,T h )ev h xi: h 

ReF(v h ,T h ) 



Vh\\DG 

Im A h (u h ,er h ;v h ,T h ) 



\vh\\DG 



sup 



lmF(v h ,T h ) 



{v h ,T h )£V h X-£ h \\vh\\DG (v h ,T h )eV h xS h \\Vh\\DG 

\F(v h ,T h )\ 

< 2 sup J-n — n 

(v h .,T h )eV h xS h \\Vh\\DG 

<Ck- x M{f,g). 



Hence (3.31) holds. 

To show (3.32), setting (v^t^) = (0, o - /,) in (3.1) and using the trace and 
Schwarz inequalities yields 

Ikftll^(n) = (Vhu h , tr H ) n + £ (iS ([\V h u h ]], [[a h ]]) e - ([[u h ]], {a h }) t 

<\\^hU h \\ 2 L 2 [n) + -\\(T h \\ 2 L2(n) +CY Yl h K 1 ( < S \\Vh'Uh\\L*(K) 

ee£l K=K e ,K' e 



+ \\Uh\\L 2 (K)^ \Wh\\L 2 {K) 



1 



Thus, 



+ C Y. J2 h K 2 (s 2 \\\/ h u h \\i 2(K) + \\u h \\ 2 L2(K) 

e£ £i h K=K C ,K' C 

r fclll=(n) £ ^ + $2 (]^ h h K 2 ))\ u h\i,h + (max /i K 2 )||w h ||^ 2(n) . 



The desired estimate ( 3.32 ) follows from combining the above inequality with (3.31 1. 



The proof is complete. 



□ 



An immediate consequence of the stability estimates is the following unique 
solvability theorem. 



Theorem 3.2. There exists a unique solution to the LDG method (3.1) for all 

k, h,S,fi> 0. 



Proof. Since problem (3.1) is equivalent to a linear system, hence, it suffices to 
show the uniqueness. But the uniqueness follows immediately from the stability 
estimates as the zero sources imply that any solution must be a trivial solution. □ 

3.2. Absolute stability of LDG method #2. In this subsection, we consider 
the LDG method #2. By adding (2.7) and (2.8) with the given numerical fluxes, 
we then recast the LDG method #2 as the following nonconforming mixed method: 
Find (lift, a%) € Vh X £/, such that 

(3.34) B h (u h ,a h ;v h ,T h ) = F(v h ,r h ) V(v h> r h ) € V h X Yi h , 
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where F is defined in (3.3 1 and 



(3.35) B h (w h , Xh\ Vh,T h ) = (x h , VhVh)n - k 2 (w h ,v h )n + lk (w h , v h ) 1 



£ ({X h }-l0[[WhM[Vh]])e 

E( 1 *<M,[[r h ]]> e -<W,W) ( 



3.2.1. A generalized inf-sup condition. The goal of this subsection is to show 



that the sesquilinear form Bh defined in (3.35) for the LDG method #2 satisfies 



another generalized inf-sup condition. To the end, we introduce the following space 
notation: 



(3.36) 



S h ■= {(w h ,Xh) € V h X £ h ; (w h ,Xh) satisfies (3.37)} 



where 



(3.37) (x h , Th)n - (VftWh, T h )n 

- E(^(WUN]) e -([W,{Ti,}) e ) -o vk,t,)6%xs,. 



We note that it is easy to check that (wh,Xh) ^ Sh implies that it satisfies (2.7) 
with uk being defined by the LDG method #2. 



Lemma 3.1. For any (wh,Xh) G Sh, there holds the following estimates: 



(3.38) \w h \ lih < \ -WXhh^n) 



2 

\L 2 (e 



s 2 
h 



ir\\[[x h }}\\l He) 



(3.39) 



\Xh - VfctUfc|| L »(n) < C(E (r-H[M]lli 2 (e) + T-||[[Xh]]|li2 (e) ^ 
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Proof. On noting that (wh,Xh) satisfies (3.37), setting = V/jW/,. in (3.37), we 

get 

\ w h\l, h = Re (X/n V h w h )n 

-Re (iHliXh]}, [[VhWh]}) e ~ ([[wh]],{VhW h }), 

+ ^2 She* \\[\Xh]]h'(e)\\VhWh\\v(K) 



ee£! 



K=K.,K' 



^he 2 Y \\il w h}}\\L2(e)\\V h W h \\ L 2 {K) 



K=K.,KL 



< \\ w h\\, h + \\\Xh\\h{U) + j^Wl,. 



+ c E(r IIMIIi2 ( e ) + r IIMIIi2 ( e 

eesf e 



Therefore, 



KlL < ^Wx h \\h m + cJ2 (^ll[K]]ll| 2(e ) + y\lix h }}\\h ( e) 



ee£! 



which gives (3.38). 



The estimate (3.39) follows from the same derivation by setting Ty L = X/j — V hWh 
in (3.37). The proof is complete. □ 

We now are ready to state a generalized inf-sup condition for the sesquilinear 
form Bfr. 

Proposition 3.2. There exits constant ci > such that for any (wh,Xh) $h 
there holds 



(3.40) 



where 



sup 



ReB h (w h) Xh> v h,T h ) 



{v h ,T h )ev h xs h \\(vh,r h )\\ 



DG 



Im B h (w h ,Xh- v h,r h ) c 2 

sup „, x,| > — \\(w h ,Xh)\\ DG , 

(v h ,T h )ev h xs h \\{vh,T h )\\ DG 72 



-k 2 + l p + s S 



(3.41) 72 : • max( , 

(3.42) HK.xOIIdg : = ( fc2 H^lli=(n)+^ 2 |l^lli=(r) 

+ \\x h \\ 2 L 2 {n) + c n \\V h w h \\ 2 L2{r) 



Proof. Since the proof follows the same lines as that of Proposition 3.1 we shall 
only highlight the main steps and point out the differences. 



Step 1: Taking the first test function. 
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We first choose the test function (vh,Th) — (wh,Xh) to get 
Bh{w h ,X h ;v h ,T h ) = B h (w h ,x h ;w h ,Xh) 

= (Xh,Xhh - k 2 {w h ,w h )n +ik (w h ,w h ) r 

E ({xJ-^[K]],[K]]) e 

'eel 

- E ( l5 (\tXh]]> llXh]})e ~ (M> {x h }), 



+ (Xh> ^hw h )n ~ (VhWh, X h )n- 
Taking the real and imaginary parts yields 

(3.43) Rc B h (w h ,Xh,Wh,Xh) = \\Xh\\h(Q) - k2 \\ w h\\l^^Q), 

(3.44) 1m B h (w h ,x h ;w h ,x h ) = k \\ w h\\ 2 L 2 {T) 



E(- 5 H[W]lli 2 (e)+^ll[W]ll! 2 (, 

'e£l 

2Im{ Xh ,V h Wh)n+ E 2Im <M> <Xj,}) e 



On noting that (u^jX/J satisfies ( |3.37[ ), setting = x,, in (3.37), we get 
IIXfcllia (n) - {V h Wh,Xhh = E ( 1(5 H[[^]]lli 2 ( e ) - ([[ w h]],{Xh}) t 



ee£! 



Taking the imaginary part yields 

Jm(x h , V h w h )n + E Im([K]HxJ) e = E 



Hence, (3.44) becomes 

(3.45) JmB h (w h ,Xh'i w h,Xh) = k \\ w h\\h(r) 

+ E(' 5 ll[[^]]lli 2 (e) + / 3 IIWII^(e: 



Step 2: Taking the second test function. 

Next, we choose another test function (vh,Th) = (ol ■ V^w^,x^), which is dif- 
ferent from the one used in the proof of Proposition |3.1[ and use the fact that 
V/jU/j = V h w h to get 

Bh(wh,Xh> Vh,r h ) = B h (w h ,x h , a ■ ^hW h ,Xh) 

= (Xh>Xh)n - k 2 (w h ,ot ■ VhWh)n + ik (w h ,a ■ V h w h ) r 
+ E i/3 UK]], [[« ■ Vhw h ]]) e - E ({Xfc}» [[« ' V h «7fc]]) e 

- E ^(M,M) e + E UKMxJ), 

+ 2ilm(x fe , VhWh). 
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Taking the real part immediately gives (vh = ol ■ V h w h) 

(3.46) ReB h (w h , Xh ;v h ,Xh) 

= \\Xh\\h(n) ~ k 2 Re(w h ,v h ) n - klm (w h ,v h ) r 

+ E ( Rc <{*J, IK]] - [M)e - Plm([[w H ]], [[v h ]]) e ). 



Step 3: Deriving an upper bound for k 2 \\wh\\ 2 L 2(^iy 

To bound (3.43) from below, we again need to get an upper bound for the term 
^ 2 1 1 w fr 1 1 i^n) on the right hand side of (3.431. 



Using the integral identity (3.10), (3.12), (3.43) and (3.46), we have 
(3.47) 2k 2 \\w h f L2{n) =2ReB h (w h , Xh ;v h ,Xh) 



(n) 



+ (d- 2) ReB h (w h ,Xh,w h ,Xh) - d \\Xh\\h 
+ 2fclm (w h ,v h ) r + k 2 (a- n n , \w h \ 2 ) r 

+ 2k 2 RcJ2 (a-n e {w h },[w h ]) e + 2ReJ2 ({X/J. M) e 
+ 2(d - 1) Re E U* J, [K]]) e + 21m £ /3 ([K]], [K]]) e 
+ d\w h \l h - E (a-n e ,|V /l w, l | 2 ) e 
-2Re E («• MVfcUto}, [V /l u; /l ]) c 
±2 Re E (« • MxJ, [V^]) e . 



We note that by (3.13) the sum of the second and third lines to the last is zero, 
and the contribution of the last line is obviously zero. These terms are purposely 
added in order to get sharper upper bounds when they are combined with the terms 
preceding them. 

We now need to bound the terms on the right-hand side of (3.47). Some of these 
have been obtained in the proof of the Proposition 3T and the others are derived 
as follows. 



(3.48) 2(d-l)Re£ <{ X J, [K]]) e 

<2(d-l)J2 h * i E \\XhhHK)\\[w h ]\\ LHe) 
<^llxJli 2( n)+CE^/3|IMIlW 
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(3.49) 2Rc £ ({ Xh }, [K]]) e - 2 Re £ (a • n e { Xh }, [V h w h ]) c 



ee£! 



ee£! 



2Re£ 



d-1 



-("•n e ){x h ■T^})Vh[u>/J 



<EE^ 5 \\X h \\v{K)\\[V h w h -Ti\\\ L * [e) 
<hx h \\l H n)+cY, ^E^II[V^-Ti]||i 2(e) . 



It follows from (3.38) and (3.39) that 



(3.50) 



2Re E ( a ' "elX/j - V^Wft}, [V/ l iy/ l ]) e 



ee £,r K=K e ,K' B 



< 



17' 



< c E(^ll[K]]lli 2(e ) + ^ll[[x/J]ll 



2 



1 

16' 



1 2 



/ 1 <5 2 

C E (j^ll[K]]lli=(e) + ^3 II Millie) 



Similar to the derivation of (3.20), we have 



(3.51) 



21m ^ /3([K]], [K]]) e = 21m £ ([w h ], [a-V h w h ]) c 

- C E 2 E ll[ U, 'i]IU 2 (e)l|V/ l W, I || L 2 (A ' ) 



K=K e ,K' 



< 1 iiv II 2 +cV^ 2 + 1 



5 2 



ll[K]]|li 2( e) + r IIMIl!2 (e) 
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Now substituting estimates ( 3.48 1~( 3.51 ) , together with (3.15), (3.16) and (3.17), 
into (3.47) we obtain 

(3.52) 2k 2 \\w h \\ 2 L2(n) < 2ReB h {w h , Xh ;v h , Xh ) 

+ (d - 2)ReB h (w h ,Xh> w h,Xh) - d \\x h \\h(n) 



+ Ck 2 \\w h f LHr) - -f W^hWhW 



2 

L2(e) 



\ w h\\mn) 



17 
1G 

1 

To 



d\\Xh\\ L 2 (Q) 



\\Xh\\ 



cdY: 



ees! 



10 

/3 2 + l 

d-1 



IIMIli» (o ) 



h 



2 

L 2 (e 



lXh\\\\L 2 (e 



-Jill?, 



lL 2 (e) 



+ ^IIX h |li»( ) + C 2 (^^IllWllll^e) + ^ll[[X/J]lli 2 (e))- 

eeel 



On noting that (3.45) provides upper bounds for terms || [[^Cft.]] Hi 2 fe) ' II t ul ' l ]lli 2 (e) 
and fc^l^D 2 2 ( r -) in terms of ImJ5^(w^, Xh ; w^, Xh)i and the jumps of the tangential 
derivatives ||[V/iW;/s • T"i]|| 2 2( e ) in (3.52) can be bounded by the inverse inequality 
(13.241), we get 



(3.53) 2fc 2 |K|| 2 2(0) + ^-\\w h \\ 2 L2{T) + ^ ]T H V ^I 



ll 2 (e) 



< 2 Re B h (w h ,Xh-, v h,Xh) + (d-2) ReB h (w h7 Xh, w h,Xh 

fc 2 ,, „o /l d 



M 3 lmB h (w h ,Xh,Wh,Xh) + '-i^\\wh\\ 2 L 2 {n y 



1G 



IXhllz,2(n) 



< 2ReB h (w h , Xh ;v h ,x h ) + (d- 1) Re B h (w h , x h \ ™h, Xh) 

3k 2 1 
+ M 3 lmB h (w h ,x h ;w h ,Xh) + -^-Ikhlli^n) ~ gll^lli 2 ^)' 



where 
(3.54) 



/fc 2 

M 3 — C[ k + max! - 



/3 + <5 <5 1 



eef^V /3/i e /i e 
Using the linearity of the sesquilinear form B^ , we have 

fe2 |Wli 2 (fi) + fc2 | Wll 2 (r) + IIX/Jli 2 (o) + cn\\V h w h \\ 2 L2(T) 

< 4ReB h (w h ,x h ;Vh,Xh) + 2 ( d ~ l )^ eB h(w h ,X^ w h,Xh) 

+ 2M 3 lmB h (w h ,Xh> w h,Xh) 
= Re B h (w h ,Xh\Wh,Xh) + 2M 3 lmB h (w h ,Xh-, w h,Xh)' 



22 



XIAOBING FENG AND YULONG XING 



where Wh = Avh + 2(d — l)wh and Vh = a- VhWh- 
Step 4- Finishing up. 

It follows from the inverse inequality (3.24) that 

k2 \\ v h\\ 2 L 2 m + k 2 \\v h \\ 2 L2(r) < Ck 2 (m&x h^ 2 \\w h \\ 2 L 2 {n) + c n \\ V h w h \\ 2 L a (r) 

Then by the definition of the norm || (wh, XiJWdg m (3-42 1 and the fact that VhVh 
VhWh, we get 



\\(.ti h ,Xh)\\ D G < 4 \\( v h,X h )\\ DG + 2{d - l)\\{w h , Xh )\ 



DG 



<C(max/i^ + fc + l) \\(w h , Xh )\\ 



DG 



Therefore, 



(3.55) 



ReB h (w h ,Xh'^h,Xh) , ImBh(wh,Xh,; w h,Xh 



> 



\{w h ,x h )\\ DG \\(wh,Xh)\\ D G 

Re B h (w h ,Xh\ Wh,Xh) 

C(maxx eTh lijl +k + l) \\{w h ,Xh)\ 



Im B h (w h ,Xh m , w h,Xh) 



DG 



> 



\\(wh,Xh)\\ 

1 ReB h (wh,Xhiu>h,Xh) + Milm B h {w h ,Xhi w h,Xh) 



DG 



2M 3 

> J \\{wh,X h )\\ DG , 

72 



\(w h ,Xh)\\ 



DG 



for some constant C2 > and 72 defined by (3.41). Hence, (3.40) holds and the 
proof is complete. □ 

3.2.2. Stability estimates. The generalized inf-sup condition proved in the last 
subsection immediately infers the following (absolute) stability and well-posedness 
theorems for the LDG method #2. 



Theorem 3.3. Let (u h ,er h ) e V h x H h solve ( |3.34[ ). Then there holds 
(3-56) ||K, < 7 h )|| DG < 72 fc- 1 M(/, 5 ). 



Proof. On noting that any solution (iih,crh) of (3.34) belongs to the set S^, the 
desired estimate (3.3) follows readily from (3.40) with (w h ,Xh) — {uinCh), (3.34) 
and (3.33). The proof is complete. □ 



Theorem 3.4. The LDG method ( 3.34 ) has a unique solution for all k, h, 5, (3 > 0. 

Since the proof of the above theorem is a verbatim copy of that of Theorem |3. 2 1 
we omit it. 



4. Error estimates 

The goal of this section is to derive error estimates for the LDG method #1 and 
#2. Following the idea of [TU], this will be done in two steps. First, we introduce an 
elliptic projection of the solution (it, er) using a corresponding coercive sesquilinear 
form of Ah (resp. Bh) and derive error bounds for the projection. We note that 
the error analysis for the elliptic projections has an independent interest in itself 
(cf. [3]). Second, we bound the error between the projection and the LDG solution 
using the stability estimates obtained in Section [3j Since the error analysis for the 
two LDG methods are similar, we shall give more details of the error analysis for 
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the LDG method #1 but shall be brief for the LDG method #2. Throughout this 
section, we let for j = 1,2 

Hj (Th) = TT H j (K), h = max h K w maxft e , /3 = /3 h~\ S = S Q h 

An KeTh ee£h 

for some positive constants /?o and So- 

4.1. Error estimates for the LDG method #1. 

4.1.1. Elliptic projection and its error estimates. For any (w,x) £ H 2 (Th) x 
-ff 1 (7/i) <i , we define the elliptic projection (tD/,,x/i) £ 14 x of (u>,x) by 

(4.1) a h {w h ,Xh\v h ,T h ) = a h (w,x;v h ,T h ) V(v h ,T h ) £ V h x ~E h , 
where 

(4.2) a h (w h ,Xh- l v h-,Th) ■ = A h (w h ,x h ; v h,r h ) + k 2 (w h ,v h ) n 

= {Xh> VhVhh + ifc (wft, «h)r 

- E ( 1<J ([[Vfcw/,]], [[Tfc]]) e - ([K]],{n}) e ) 

To derive error bounds for the above elliptic projection, we first notice that 
a h (w h , x h - Vh, -VhVh) = (VhWh,V h v h ) n + lk (wh, v h ) r 

- E ({V h w h }-ip h- 1 [[w h ]},[[v h }}) e 

+ E («J ([[Vfc^]], [[v^]]> e - ([K]], {v hUfc » e ) 

= : Ah(w h ,v h ). 
As a result, e V/, satisfies 

(4.3) = A h (w,v h ) Vv h e Vfc. 
Moreover, since 

E (^([[v^]],M) e - <W,{T fc }) e ), 

we have that %h £ S/, satisfies 
(4-4) (x h , r h ) n = {V h w h , r h )n 

+ E (iS ([N h w h - V h w}], [[r h ]]) e ([[w h - w}], {r h }) e ) 

+ {x, T h)n - (Vfciu, T ft )o Vt/» e E h . 
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Lemma 4.1. For any w,v £ H 2 (Tii), there exists a k- and h-independent constant 
C such that 



(4.5) 



\A h (w,v)\ < cilHllli^HHlli^. 



Moreover, for any e £ (0, 1), there exists a constant c e > such that 



(4.6) ReA h (v h ,v h ) + {l-e + c £ )lmA h (v h ,v h ) > (1 - e)\\v h \\l h , 



whe 



(4.7) 



(4.8) 



HHkh : = {\\^hw\\ 2 L 2 [n) + k\\w\\ 2 L2{T) 

+ E(' 9 IINI| 2 L 2 (e)+'5|l[[VH]l| 2 L 2 ( e ) 



ee£, 



\ w \\\l 



, h -= (iH^+E^IKv^-ne}!! 



2 

L 2 « 



Since the proof of the above lemma is elementary, we omit it. We now recall the 
following stability estimate for u (cf. [71 ITU]): 

IHI^cnj^^ + W/.fl), 

which is needed to prove the next lemma and will be used several times in the rest 
of this section. 



Proposition 4.1. Let u € H 2 (£l) be the solution to problem (1.1) -(1.2) and 
a = Vu. Let (uf l ,d'f l ) £ Vh x denote the elliptic projection o/(u, <r) defined by 



(4.1). Then there hold the following error estimates: 



(4.9) 

(4.10) 

(4.11) 



\u-u h \\i,h + k^\\u-u h \\ L 2 (r) < (l + kh)ikh, 
\\u - u h \\ L 2 {n) < (1 + kh)kh 2 , 
Ik - a- h \\ L 2 {n) < (1 + kh)^kh. 



Proof. Since the proof of (4.9) and (4.10) is essentially same as that of [TOl Lemma 
5.2], we omit it to save the space and refer the reader to [TU] for the details. 

To show (4.11), on noting that (4.4) and the identity (a,Th) = (VhU,Th) imply 



(4.12) ((T-a h ,T h ) n = (V h u - V h fih, T h )n - ^ fah{[\V h (u-u h )]][[T h ]]) 



([[u-U h ]],{T h }) e ) VT h £-Z h 
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For any Xh S we set t/, = ^ — <t/j. Then by (4.12), the trace inequality, 
Schwarz inequality we get 

Ik - <7h\\L 2 (n) = { <T — <Th,cr - Xh)ct + {<r — <Jh,T~h)n 
= (a - cf h ,cr - Xh)n + (V h (it - u h ),T h ) n 

- J2 fati ([[V h (u - u h )]}, [[r h ]]) e - ([[u - u h ]}, {th}), 
eee' 



< || tr - CTh|| L 2(n)||er - XhW^(Q) + \\^h(u - u h )\\ L 2( n) \\T h \\ L 2(n) 

+ C(J2 Ml t^C" - fifc)]]lli»(e) + [« - fi/JI&Ce) 



■ £ 



2 

L 2 (e 



< 



jlk - CTh||i2 (n) + ^llT/Jl^fo) + ||<x - X/,llia(n) + c \\ u - Uh\\i,h 



' - II - <T,,ll7 . on + 7||CT-X ft ||| 2(n) +C||u- U h \\i, h . 



I CT - °"/»ll£»(n) + -| 



Hence, it follows from the above inequality, (4.9 ), and the polynomial approximation 
theory (cf. [2]) that 

Ik - ^|| L 2 (n) < C||w-%||i )ft + 2 inf \\(T - Xh\\^{n) 

Xh^h 

< {l + kh^hk+ik + k-^h 

< (l + kh)?hk, 



which gives (4.11). The proof is complete. 



□ 



4.1.2. Global error estimates for the LDG method #1. In the preceding 
subsection we have derived the error bounds for (u — Uh, cr — cr^). By the decom- 
position u - u h = (u - u h ) + (u h - Uh) and cr - cr h = (cr - cr h ) + {& h - cr h ) and 
the triangle inequality, it suffices to get error bounds for [uu — Uhi&h — Vh)- We 
shall accomplish this task by exploiting the linearity of the Helmholtz equation and 
using the stability estimate for the LDG method #1 obtained in Section |3~T) 
First, on noting that (u, cr) satisfies 



(4.13) 



A h (u,a-,Vh,T h ) = F(v h ,Th) V(w /i ,t /i ) G V h x E h . 



Subtracting (3.1) from ( 4.13| ) yields the following error equation (or Galerkin or- 
thogonality) : 

(4.14) A h (u-u h ,er - a- h ;v h ,Th) =0 V(v h ,T h ) e V h x T, h . 

Next, to proceed we introduce the notation 

u-u h = e h + q h , e h :=u-u h , q h :=u h -u h , 

cr- cr h =ip h + <j> h , tPh'-= CT - °7i) <t>h '■= °h - cr h . 
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Then by (4.14) and the definitions of the sesquilinear form at and the elliptic 
projection we have 

(4.15) A h (q h) ct> h ;v h ,T h ) = -A h (e h ,tp h ;v h ,T h ) 

= -a h (e h ,iJj h ;v h ,Th) + k 2 (e hl v h ) 

= k 2 (e h ,v h ), V(v h ,T h ) e V h x H h . 

The above equation implies that (qh,(f>h) € Vh x is the LDG solution to the 
Helmholtz problem with source terms / = k 2 eh and g = 0. Then an application of 
the stability estimates of Theorem |3.1| immediately yields the following lemma. 



Proposition 4.2. There hold the following estimates for 

(4.16) II^Hdg <7i(l + ^)fc 2 ^ 2 , 

(4.17) \\4>h\\v{tt) £ 7i(l + W + S kh 2 )(l + kh)kh. 



Combining Proposition 4.1 and 4.2 using the triangle inequality and the stan- 
dard duality argument give the following error estimates for (v,h, &h)- 

Theorem 4.1. Let u € H 2 (fl) be the solution to problem — ( 1.2 ) and a := Vu, 
and (uh,(7h) be the solution to problem (3.1). Then there hold the following error 
estimates for (uh, gh)'- 



(4.18) \\u - u h \\ hh + ki\\u - u h \\ L 2 (r) < ((l + kh)i +^(1 + kh)kh)kh, 

(4.19) ||u - u h \\ L 2 (n) < (1 + 7l )(l + kh)kh 2 , 

(4.20) \\tr - o-h|| za( n) < ((1 + kh)i + 7l (l + kh)(l + kh + 5 a kh 2 ))kh. 

4.2. Error estimates for the LDG method #2. The error analysis for the 
LDG method #2 essentially follows the same lines as that for the LDG method #1 
given in the previous subsection. However, there are three main differences which 
we now explain. First, the sesquilinear form ah needs to be replaced by another 
sesquilinear form bh in the definition of the elliptic projection (4.1), where bh is 
defined by 



(4-21) 



bh(wh,Xh> v h,Th) ■ = B h (w h ,Xh'i v h,T h ) + k 2 (w hl v h )n 
= (Xhi ^hVh)n + ifc (w h , v h ) r 

-E«^}-i/3[K]],[K]]> e 



E 

eeel 



(^<[[xJ],[[r fe ]]) e -(W,{r,}>, 



(Xh> T h)n - (V h w hl T h ) 



Second, due to strong coupling between Uh and a hi the error estimates for the 
new elliptic projection (uh,&h) niust be derived differently. To the end, we need 
the following lemma, which replaces Lemma |4.1| 



Lemma 4.2. Let /3 — fioh^ 1 and 5 — 5$h for some positive constants (3q and Sq. 

(i) There exists an h- and k-independent constant C3 > such that the sesquilin- 
ear form bh satisfies the following generalized inf-sup condition: for any fixed 
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(4.22) sup ^^K.Xhi^.Th) 

(v h ,T h )eV h xS h \\\{vh,r h )\\\ D G 



+ sup 

(«h,T h )ev r h xS J 



> C 3 \\\(w h ,Xh)\\\DG- 

Vh,Th)\\\DG 

(ii) There exists an h- and k -independent constant C > such that for any 
(w, X ), («,r) G i? 2 (T ft ) x H 1 (T h ) d , there holds 

b h (w,x;v,T)\ < C\\\{w,x) 



(4.23) 

where 

(4.24) 
(4.25) 



{w,x) 



Il 2 (o) 



Idg 



/^IIMI' 2 



L 2 (e 



The proof of (i) is based on evaluating the first quotient on the left-hand side of 
(4.22) at (vh,Th) = ((l + Ci)iVh,CiXh~^h w h) and evaluating the second quotient 
at (vh, Th) — {C2Wh,C2Xh) f° r some sufficiently large positive constants C\ and C%. 
The proof of (ii) is a straightforward application of Schwarz and trace inequalities. 
We skip the rest of the derivation to save space. 

The above generalized inf-sup condition, the boundedness of the sesquilinear 
form bh, and the duality argument (cf. [5j) readily infer the following error estimates 
for the new elliptic projection (u^, er/j). We omit the proof since it is standard. 



Proposition 4.3. Under the assumptions of Proposition J^.l . there hold the fol- 
lowing estimates: 

(4.26) \\u - u h \\i,h + \\<r - 6Vi|U=>(fi) < kh, 

(4.27) \\u-u h \\ L 2 [n) <k 2 h 2 . 

The third difference is that the new error function (qh, (p^) now satisfies 

(4.28) B h (q h ,(p h ;v h ,T h ) = k 2 (e h ,v h ) V(v h , r h ) G V h x S h . 



As a result, by Theorem 3.3 and (3.38) we get 

\q h \l,h + \\{q h ,cj) h )\\ DG < 72(1 + kh)k 2 h 2 , 



(4.29) 



which replaces estimates (4.16) and (4.17). 



After having established Proposition 4.3 and (4.29), once again, by the triangle 
inequality we arrive at the following error estimates for the solution (uh,crh) to the 
LDG method #2. 



Theorem 4.2. LetueH 2 (tt) be the solution to problem (1.1) -(1.2) andcr:=Vu, 
and (uh,(Th) be the solution to problem (3.34). Then there hold the following error 
estimates for (uh, (Th)'- 

(4.30) \\u - u h \\i,h + kS \\u - u h ||i2(r) 

+ ||<7 - <T h \\ L 2 (n) < (l + 72(1 + kh)kh))kh, 

(4.31) ||u - u h \\ L 2 (n) < (1 + 72 (1 + kh))k 2 h 2 . 
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Remark 4.1. (4.29) shows that <p h := er^ — <Th has an optimal order (in h) error 
bound for the LDG method #2, while (4.17) shows that <p h only has a sub-optimal 
order error bound for the LDG method #1. We believe that this is the main reason 
why in practice the LDG method #2 gives a better approximation to the flux variable 
<r than the LDG method #1 does although both methods have the same asymptotic 
rate of convergence in h. 

5. Numerical experiments 

In this section we shall provide some numerical results of the two proposed LDG 
methods. Our tests are done for the following 2-d Helmholtz problem: 

sin(fcr) 



(5.1) 
(5.2) 



-Au - k z u = f := 
du 

lku = < 



in Q, 
on := dfl. 



Here fi is the unit square [—0.5, 0.5] x [—0.5, 0.5], and g is chosen so that the exact 
solution is given by 

cos(fcr) cosfc + isinfc 



(5.3) 



-J (kr) 



k k(j (k)+iJx(k))' 

in polar coordinates, where J v (z) are Bessel functions of the first kind. 

Assume T\/ m be the regular triangulation that consists of 2m 2 right-angled 
equicrural triangles of size h = 1/m, for any positive integer m. See Figure [l] 
for the sample triangulation T1/4 and 7i/io- 




Figure 1. The computational domain and sample meshes. Left: 
7i/4 that consists of right-angled equicrural triangles of size h = 4; 
Right: T[/io with h — ^. 



5.1. Sensitivity with respect to the parameters 8 and (3. In this subsection, 
we examine the sensitivity of the error of the LDG solutions in H^-seminorm with 
respect to the parameters 5 and (3. 

The LDG method #1 is considered first. We start by fixing 5 — 0.1h e and 
testing the sensitivity in the parameter f3. With two wave numbers k = 5 and 
50, we compute the solutions of the LDG method #1 with different values of j3: 
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O.OOlh^ 1 , O.Qlh^ 1 , /ijT 1 and 1. The relative errors, denned by the errors in the 
-ff^-seminorm divided by the exact solution in the 77 1 -seminorm, are shown in the 
left graph of Figure [2] We observe that the relative errors have similar behaviors 
and decay as mesh size h becomes smaller. This shows that the errors are not 
sensitive to the parameter 0. Next, we fix f3 = 0.001/iJ 1 , and repeat the test with 
different S. The right graph of Figure [2] shows the relative errors with parameters 
8 = 0.001/i e , 0.lh e , 10h e and 0.1, and wave numbers k — 5 and 50. We observe that 
the errors have similar behaviors for small values 6 = 0.001/i e and 0.1h e . Larger 6 
results in larger error. 

The sensitivity tests of the LDG method #2 are shown in Figure [3j similar 
behaviors are also observed. 







k=50 


Solid: p=0.001/h 
_ Dashed: p=0.01/h 
Dotted: [S=1 
Dashdot: |B=1/h 







« k=S 




Solid: 6=0.1 *h 
Dotted: 5=0.001*: 
Dashed: 5=0.1 
Dashdot: 5=1 0*h 



Figure 2. Relative error in the if 1 -seminorm of the LDG method 
#1 with different parameters for two wave numbers k = 5 and 50. 
Left: S = 0.1h e is fixed, /3 = O.OOl/i^ 1 , OMh- 1 , hj 1 and 1; Right: 
/3 = O.OOl/i^ 1 is fixed, 6 = 0.0Ql/i e , 0.1h e , 10h e and 0.1. 



5.2. Errors of the LDG solutions. In this subsection, we fix the parameters and 
investigate the changes of the numerical errors as functions of the mesh size. 

We start from the LDG method #1. As suggested by the sensitivity tests in the 
previous subsection, we pick 



(5.4) 



0.lh e 



13 = 0.001/1" 



The relative error of the LDG method, and the finite element interpolation are 
shown in the left graph of Figure |4j with four different wave numbers k = 5, 10, 
50 and 100. The relative error of the LDG solution stays around 100% before a 
critical mesh size is reached, then decays at a rate greater than —1 in the log-log 
scale but converges as fast as the finite element interpolation (with slope —1) for 
small h. The critical mesh size decreases as k increases. 

The right graph of Figure [4] contains the relative error when we fix kh = 1 and 
hk = 0.5. It indicates that unlike the error of the finite element interpolation the 
error of the LDG is not controlled by the magnitude of kh, which suggests that there 
is a pollution contribution in the total error. The left graph of Figure [5] contains 
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k=50 


Solid: p=0.001/h 
_ Dashed: p=0.01/h \. 
Doited: p=1 
Dashdot: |B=1/h 




k=5 




Solid: 6=0.1 *h 
Dotted: 5=0.001*h 
Dashed: 5=0.1 
Dashdot: 5=1 0*h 



Figure 3. Relative error in the i? 1 -seminorm of the LDG method 
7^2 with different parameters for two wave numbers k = 5 and 50. 
Left: 5 = 0.1h e is fixed, /3 = 0.001/i~\ O.Ol/i^ 1 , li; 1 and 1; Right: 
(3 = O.OOl^ 1 is fixed, <5 = 0.001/i e , 0.1/i e , 10/i e and 0.1. 



the relative error of the LDG method with the mesh size satisfying k 3 h 2 = 1 for 
different values of h. The error does not increase with respect to k. 

The LDG method #2 has also been tested using the same parameters in (5.4). 
Similar behaviors have been observed as shown in Figure [5] and [6j 







k=loo 




"V^ k=50 




^S5^N k=10 





^ k=5 











kh=1 




kh=0.5 


: 1 






Figure 4. Left: relative error of the LDG method #1 (solid 
line) and the finite element interpolation (dotted line) in the H 1 - 
scminorm for k = 5, 10, 50 and 100; Right: relative error of the 
LDG method #1 in the i? 1 -seminorm for k = 1, • • • , 200, kh = 1 
and kh — 0.5. 

At the end, we look closely at the situation with a large relative error when 
kh > 1. The LDG method #1 solution with parameters <5 = 0.1h e , (3 — O.OOl/i" 1 , 
k = 100 and h — 1/45 has a large relative error of size 0.9392. The surface plots 
of the finite element interpolation and the LDG solution are given in Figure [7J It 
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0.5 I 1 10 1 



0.4 




FIGURE 5. Left: relative error of the LDG method #1 in the H 1 - 
seminorm with k 3 h 2 — 1; Right: relative error of the LDG method 
#2 (solid line) and the finite element interpolation (dotted line) in 
the i? 1 -seminorm for k = 5, 10, 50 and 100. 




shows that the LDG solution has the correct shape/phase although its amplitude 
is smaller. 

5.3. Comparison between the two LDG methods. Two different LDG meth- 
ods are proposed in this paper. The first one is derived following the IPDG method 
proposed in [10] . and the second one has a more standard numerical flux formula- 
tion and is supposed to have a better approximation for the vector/flux variable. 
In this subsection, we provide a comparison between these two methods, in terms 
of the error and computational cost. 

We start by revisiting the test examples of Subsection |5.1[ Instead of computing 
the relative error of Uh in the i? 1 -seminorm, we compute the relative error of cr^ in 
the L 2 -norm for the LDG method #2. The numerical results are presented in Figure 
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Figure 7. Left: surface plots of the finite element interpolation 
(left) and the LDG method #1 solution (right) with parameters 
S = 0.1h e , /3 = O.OOl/i" 1 , k = 100 and h = 1/45. 

[8j which show that although the solution is still not sensitive to the parameter /?, 
better approximation to cr is achieved for larger 5. It confirms our prediction that 
the LDG method #2 gives a better approximation for the vector/flux variable. 
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Figure 8. Relative error of crt in L 2 -norm of the LDG method 
7^2 with different parameters for two wave numbers k = 5 and 50. 
Left: S = 0.1h e is fixed, (3 = 0.001/j- 1 , O.Ol/i,: 1 , h- 1 and 1; Right: 
13 = O.OOl^ 1 is fixed, <5 = 0.001/i e , 0.1h e , 10h e and 0.1. 

Table [l] provides a detailed comparison of these two methods for different mesh 
sizes h, with the parameters S — 0.lh e , (3 — 0.001/iJ 1 and k — 10. It shows that 
the computational cost of the LDG method #2 is about twice larger than that of 
the LDG method #1. Also, as expected, the error of the vector/flux variable of the 
LDG method #2 is smaller than that of the LDG method #1, and both methods 
demonstrate a first order rate of convergence. 



5.4. Comparison between LDG and finite element solutions. We have shown 
the performance and comparison of the two LDG methods in previous subsections. 
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Table 1. Comparison of the two LDG methods with parameters 
S = 0.1h e and /3 = 0.001/iJ 1 . 





1/h 


\u-u h \ H i 


order 


Ik - (r h \\ L 2 


order 


CPU time (s) 


LDG #1 


5 


4.1059E-01 




5.4715E-01 




0.0641 


10 


1.6915E-01 


1.2794 


2.4712E-01 


1.1467 


0.2381 


20 


7.6089E-02 


1.1525 


1.1804E-01 


1.0659 


0.9671 


40 


3.6648E-02 


1.0539 


5.7114E-02 


1.0474 


3.9380 


80 


1.8151E-02 


1.0137 


2.8319E-02 


1.0121 


15.8194 


160 


9.0379E-03 


1.0060 


1.4004E-02 


1.0159 


69.1861 


LDG #2 


5 


2.4711E-01 




2.2184E-01 




0.1057 


10 


1.4040E-01 


0.8156 


7.6775E-02 


1.5308 


0.4368 


20 


6.6992E-02 


1.0675 


3.3630E-02 


1.1909 


1.8356 


40 


3.2693E-02 


1.0350 


1.5710E-02 


1.0981 


7.9357 


80 


1.6165E-02 


1.0161 


7.7418E-03 


1.0209 


34.5688 


160 


8.0949E-03 


0.9978 


3.9127E-03 


0.9845 


157.6018 




Figure 9. The traces of the LDG #1 solution (left) and the finite 
element solution (right) in the xz-plane, for k = 100 and h = 
1/50 (top), 1/120 (middle) and 1/200 (bottom), respectively. The 
dotted lines are the traces of the exact solution. 



In this subsection, we provide a brief comparison between the LDG solution and 
the Pi conforming finite element solution. 
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Wc consider the Hclmholtz problem (5.1 1-(5.2 1 with wave number k = 100. With 
mesh size h = 1/50, 1/120 and 1/200, we plot the traces of the LDG method #1 
solution with parameters (5.4) in xz-plane in the left column of Figure[9j The exact 
solution is also provided as a reference. The traces of the finite element solution 
are shown in the right column of Figure [9] It is clear that the LDG method #1 
has a better approximation to the exact solution. Larger phase error in the finite 
element solution is observed in all three cases. Also, the LDG solution has a better 
approximation for the amplitude of the exact solution. 
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